#include <iostream>
#include "IsoBandPolygons.h"
#include "OCMain/Timer.h"

int computeCenterAverage(double x0, double x1, double x2, double x3, double minV, double maxV)
{
    double avg = (x0 + x1 + x2 + x3) / 4.0;
    if (avg < minV)
        return 0;
    if (avg > maxV)
        return 2;
    return 1;
}

class ShapeCoordinates
{
public:
    static void square() {}
    static void triangle_bl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate(x0, x1, minV, maxV);
        double leftbottom = interpolate(x0, x3, minV, maxV);
        cell.edges.data["lb"] = {{
                                     0,
                                     leftbottom,
                                     bottomleft,
                                     0,
                                 },
                                 0,
                                 -1,
                                 "tl"};
    }
    static void triangle_br(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);
        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                1,
                rightbottom,
            },
            1,
            0,
            "lb"};
    }
    static void triangle_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double righttop = interpolate(x1, x2, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);
        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                topright,
                1,
            },
            0,
            1,
            "br"};
    }
    static void triangle_tl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate(x3, x2, minV, maxV);
        double lefttop = interpolate(x0, x3, minV, maxV);
        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                0,
                lefttop,
            },
            -1,
            0,
            "rt"};
    }
    static void tetragon_bl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                leftbottom,
            },
            -1,
            0,
            "rb"};

        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                bottomright,
                0,
            },
            0,
            -1,
            "tr"};
    }
    static void tetragon_br(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                1,
                righttop,
            },
            1,
            0,
            "lt"};

        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomright,
                0,
            },
            0,
            -1,
            "tr"};
    }
    static void tetragon_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);
        cell.edges.data["rb"] = {
            {1, rightbottom,
             topleft, 1},
            0,
            1,
            "bl",
        };

        cell.edges.data["tr"] = {
            {topright, 1,
             1, righttop},
            1,
            0,
            "lt",
        };
    }
    static void tetragon_tl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);

        cell.edges.data["tr"] = {
            {topright, 1,
             0, leftbottom},
            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {0, lefttop,
             topleft, 1},
            0,
            1,
            "bl",
        };
    }
    static void tetragon_b(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);
        cell.edges.data["lb"] = {
            {0, leftbottom,
             1, rightbottom},
            1,
            0,
            "lb"};
    }
    static void tetragon_r(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);
        cell.edges.data["br"] = {
            {bottomright, 0, topright, 1},
            0,
            1,
            "br"};
    }
    static void tetragon_t(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double righttop = interpolate(x1, x2, minV, maxV);
        double lefttop = interpolate(x0, x3, minV, maxV);
        cell.edges.data["rt"] = {
            {1, righttop,
             0, lefttop},
            -1,
            0,
            "rt",
        };
    }
    static void tetragon_l(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate(x3, x2, minV, maxV);
        double bottomleft = interpolate(x0, x1, minV, maxV);

        cell.edges.data["tl"] = {
            {topleft, 1,
             bottomleft, 0},
            0,
            -1,
            "tl",
        };
    }
    static void tetragon_lr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                1,
                righttop,
            },

            1,
            0,
            "lt",
        };

        cell.edges.data["rb"] = {
            {1, rightbottom,
             0, leftbottom},
            -1,
            0,
            "rb",
        };
    }
    static void tetragon_tb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double bottomleft = interpolate_a(x0, x1, minV, maxV);

        cell.edges.data["tr"] = {
            {
                topright,
                1,
                bottomright,
                0,
            },

            0,
            -1,
            "tr",
        };

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                topleft,
                1,
            },
            0,
            1,
            "bl",
        };
    }
    static void pentagon_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate(x3, x2, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);

        cell.edges.data["tl"] = {
            {topleft, 1, 1,
             rightbottom},

            1,
            0,
            "lb",
        };
    }
    static void pentagon_tl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                topright,
                1,
            },
            0,
            1,
            "br",
        };
    }
    static void pentagon_bl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double lefttop = interpolate(x0, x3, minV, maxV);
        double bottomright = interpolate(x0, x1, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                0,
                lefttop,
            },
            -1,
            0,
            "rt",
        };
    }
    static void pentagon_br(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate(x0, x1, minV, maxV);
        double righttop = interpolate(x1, x2, minV, maxV);

        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                bottomleft,
                0,
            },
            0,
            -1,
            "tl",
        };
    }
    static void pentagon_tr_rl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double lefttop = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                1,
                righttop,
            },
            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                0,
                lefttop,
            },
            -1,
            0,
            "rt",
        };
    }
    static void pentagon_rb_bt(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double righttop = interpolate(x1, x2, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);

        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                bottomright,
                0,
            },

            0,
            -1,
            "tr",
        };
        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                topright,
                1,
            },

            0,
            1,
            "br",
        };
    }
    static void pentagon_bl_lr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                0,
                leftbottom,
            },
            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                1,
                rightbottom,
            },
            1,
            0,
            "lb",
        };
    }
    static void pentagon_lt_tb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double bottomleft = interpolate(x0, x1, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                topleft,
                1,
            },

            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                bottomleft,
                0,
            },

            0,
            -1,
            "tl",
        };
    }
    static void pentagon_bl_tb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double lefttop = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate(x3, x2, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double bottomleft = interpolate_a(x0, x1, minV, maxV);

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                lefttop,
            },

            -1,
            0,
            "rt",
        };
        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                bottomright,
                0,
            },

            0,
            -1,
            "tr",
        };
    }
    static void pentagon_lt_rl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate(x1, x3, minV, maxV);

        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                topright,
                1,
            },
            0,
            1,
            "br",
        };
        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                0,
                leftbottom,
            },

            -1,
            0,
            "rb",
        };
    }
    static void pentagon_tr_bt(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);
        double bottomright = interpolate(x0, x1, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                topleft,
                1,
            },

            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                1,
                rightbottom,
            },

            1,
            0,
            "lb",
        };
    }
    static void pentagon_rb_lr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);
        double bottomleft = interpolate(x0, x1, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                1,
                righttop,
            },
            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomleft,
                0,
            },
            0,
            -1,
            "tl",
        };
    }
    static void hexagon_lt_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                topleft,
                1,
            },

            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                1,
                rightbottom,
            },
            1,
            0,
            "lb",
        };
    }
    static void hexagon_bl_lt(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                0,
                leftbottom,
            },

            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                topright,
                1,
            },
            0,
            1,
            "br",
        };
    }
    static void hexagon_bl_rb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double lefttop = interpolate(x0, x3, minV, maxV);
        double righttop = interpolate(x1, x2, minV, maxV);

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                lefttop,
            },
            -1,
            0,
            "rt",
        };
        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                bottomright,
                0,
            },
            0,
            -1,
            "tr",
        };
    }
    static void hexagon_tr_rb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate(x0, x1, minV, maxV);
        double topleft = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                1,
                righttop,
            },

            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomleft,
                0,
            },
            0,
            -1,
            "tl",
        };
    }
    static void hexagon_lt_rb(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate(x1, x2, minV, maxV);
        double bottomleft = interpolate(x0, x1, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                topright,
                1,
            },
            0,
            1,
            "br",
        };
        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                bottomleft,
                0,
            },
            0,
            -1,
            "tl",
        };
    }
    static void hexagon_bl_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double lefttop = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate(x3, x2, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                0,
                lefttop,
            },
            -1,
            0,
            "rt",
        };
        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                1,
                rightbottom,
            },

            1,
            0,
            "lb",
        };
    }
    static void heptagon_tr(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double topright = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate(x1, x2, minV, maxV);

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                leftbottom,
            },

            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                topright,
                1,
            },
            0,
            1,
            "br",
        };
        cell.edges.data["rt"] = {
            {
                1,
                righttop,
                bottomright,
                0,
            },
            0,
            -1,
            "tr",
        };
    }
    static void heptagon_bl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate(x0, x1, minV, maxV);
        double leftbottom = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["lb"] = {
            {
                0,
                leftbottom,
                topleft,
                1,
            },

            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                1,
                righttop,
            },
            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomleft,
                0,
            },
            0,
            -1,
            "tl",
        };
    }
    static void heptagon_tl(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double lefttop = interpolate(x0, x3, minV, maxV);
        double topleft = interpolate(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                lefttop,
            },
            -1,
            0,
            "rt",
        };
        cell.edges.data["tl"] = {
            {
                topleft,
                1,
                1,
                righttop,
            },
            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomright,
                0,
            },
            0,
            -1,
            "tr",
        };
    }
    static void heptagon_br(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomright = interpolate(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double rightbottom = interpolate(x1, x2, minV, maxV);

        cell.edges.data["br"] = {
            {
                bottomright,
                0,
                0,
                leftbottom,
            },

            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {0, lefttop, topleft, 1},
            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                1,
                rightbottom,
            },

            1,
            0,
            "lb",
        };
    }
    static void octagon(Cell &cell, double x0, double x1, double x2, double x3, double minV, double maxV)
    {
        double bottomleft = interpolate_a(x0, x1, minV, maxV);
        double bottomright = interpolate_b(x0, x1, minV, maxV);
        double leftbottom = interpolate_a(x0, x3, minV, maxV);
        double lefttop = interpolate_b(x0, x3, minV, maxV);
        double topleft = interpolate_a(x3, x2, minV, maxV);
        double topright = interpolate_b(x3, x2, minV, maxV);
        double righttop = interpolate_b(x1, x2, minV, maxV);
        double rightbottom = interpolate_a(x1, x2, minV, maxV);

        cell.edges.data["bl"] = {
            {
                bottomleft,
                0,
                0,
                leftbottom,
            },

            -1,
            0,
            "rb",
        };
        cell.edges.data["lt"] = {
            {
                0,
                lefttop,
                topleft,
                1,
            },

            0,
            1,
            "bl",
        };
        cell.edges.data["tr"] = {
            {
                topright,
                1,
                1,
                righttop,
            },

            1,
            0,
            "lt",
        };
        cell.edges.data["rb"] = {
            {
                1,
                rightbottom,
                bottomright,
                0,
            },
            0,
            -1,
            "tr",
        };
    }
};

class IsoBand
{
public:
    static IsoBand *createInstance()
    {
        return new IsoBand();
    }

    static void releaseInstance(IsoBand *ptr)
    {
        delete ptr;
    }

    Cell prepareCell(std::vector<double> &data, int width, int x, int y, double minV, double maxV)
    {
        Cell cell;
        cell.x = x;
        cell.y = y;

        /*  compose the 4-trit corner representation */
        int cval = 0;
        double x0 = data[y * width + x];
        double x1 = data[y * width + x + 1];
        double x2 = data[(y + 1) * width + x + 1];
        double x3 = data[(y + 1) * width + x];

        cell.x0 = x0;
        cell.x1 = x1;
        cell.x2 = x2;
        cell.x3 = x3;

        if (std::isnan(x0) || std::isnan(x1) || std::isnan(x2) || std::isnan(x3))
        {
            return cell;
        }

        cval |= x3 < minV ? 0 : x3 > maxV ? 128
                                          : 64;
        cval |= x2 < minV ? 0 : x2 > maxV ? 32
                                          : 16;
        cval |= x1 < minV ? 0 : x1 > maxV ? 8
                                          : 4;
        cval |= x0 < minV ? 0 : x0 > maxV ? 2
                                          : 1;

        /* make sure cval is a number */
        cell.cval = cval;

        /*
         * cell center average trit for ambiguous cases, where
         * 0 ... below iso band
         * 1 ... within iso band
         * 2 ... above isoband
         */
        int center_avg = 0;

        /*
         * Compute interpolated intersections of the polygon(s)
         * with the cell borders and (i) add edges for polygon
         * trace-back, or (ii) a list of small closed polygons
         * according to look-up table
         */
        switch (cval)
        {
        case 85 /* 1111 */:
            ShapeCoordinates::square();
        /* fall through */
        case 0: /* 0000 */
        /* fall through */
        case 170 /* 2222 */:
            break;

            /* single triangle cases */

        case 169 /* 2221 */:
            ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 166 /* 2212 */:
            ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 154 /* 2122 */:
            ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 106 /* 1222 */:
            ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 1 /* 0001 */:
            ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 4 /* 0010 */:
            ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 16 /* 0100 */:
            ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 64 /* 1000 */:
            ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

            /* single trapezoid cases */

        case 168 /* 2220 */:
            ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 162 /* 2202 */:
            ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 138 /* 2022 */:
            ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 42 /* 0222 */:
            ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 2 /* 0002 */:
            ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 8 /* 0020 */:
            ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 32 /* 0200 */:
            ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 128 /* 2000 */:
            ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

            /* single rectangle cases */

        case 5 /* 0011 */:
            ShapeCoordinates::tetragon_b(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 20 /* 0110 */:
            ShapeCoordinates::tetragon_r(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 80 /* 1100 */:
            ShapeCoordinates::tetragon_t(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 65 /* 1001 */:
            ShapeCoordinates::tetragon_l(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 165 /* 2211 */:
            ShapeCoordinates::tetragon_b(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 150 /* 2112 */:
            ShapeCoordinates::tetragon_r(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 90 /* 1122 */:
            ShapeCoordinates::tetragon_t(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 105 /* 1221 */:
            ShapeCoordinates::tetragon_l(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 160 /* 2200 */:
            ShapeCoordinates::tetragon_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 130 /* 2002 */:
            ShapeCoordinates::tetragon_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 10 /* 0022 */:
            ShapeCoordinates::tetragon_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 40 /* 0220 */:
            ShapeCoordinates::tetragon_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

            /* single pentagon cases */

        case 101 /* 1211 */:
            ShapeCoordinates::pentagon_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 149 /* 2111 */:
            ShapeCoordinates::pentagon_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 86 /* 1112 */:
            ShapeCoordinates::pentagon_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 89 /* 1121 */:
            ShapeCoordinates::pentagon_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 69 /* 1011 */:
            ShapeCoordinates::pentagon_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 21 /* 0111 */:
            ShapeCoordinates::pentagon_tl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 84 /* 1110 */:
            ShapeCoordinates::pentagon_bl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 81 /* 1101 */:
            ShapeCoordinates::pentagon_br(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 96 /* 1200 */:
            ShapeCoordinates::pentagon_tr_rl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 24 /* 0120 */:
            ShapeCoordinates::pentagon_rb_bt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 6 /* 0012 */:
            ShapeCoordinates::pentagon_bl_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 129 /* 2001 */:
            ShapeCoordinates::pentagon_lt_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 74 /* 1022 */:
            ShapeCoordinates::pentagon_tr_rl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 146 /* 2102 */:
            ShapeCoordinates::pentagon_rb_bt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 164 /* 2210 */:
            ShapeCoordinates::pentagon_bl_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 41 /* 0221 */:
            ShapeCoordinates::pentagon_lt_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 66 /* 1002 */:
            ShapeCoordinates::pentagon_bl_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 144 /* 2100 */:
            ShapeCoordinates::pentagon_lt_rl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 36 /* 0210 */:
            ShapeCoordinates::pentagon_tr_bt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 9 /* 0021 */:
            ShapeCoordinates::pentagon_rb_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 104 /* 1220 */:
            ShapeCoordinates::pentagon_bl_tb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 26 /* 0122 */:
            ShapeCoordinates::pentagon_lt_rl(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 134 /* 2012 */:
            ShapeCoordinates::pentagon_tr_bt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 161 /* 2201 */:
            ShapeCoordinates::pentagon_rb_lr(cell, x0, x1, x2, x3, minV, maxV);
            break;

            /* single hexagon cases */

        case 37 /* 0211 */:
            ShapeCoordinates::hexagon_lt_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 148 /* 2110 */:
            ShapeCoordinates::hexagon_bl_lt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 82 /* 1102 */:
            ShapeCoordinates::hexagon_bl_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 73 /* 1021 */:
            ShapeCoordinates::hexagon_tr_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 133 /* 2011 */:
            ShapeCoordinates::hexagon_lt_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 22 /* 0112 */:
            ShapeCoordinates::hexagon_bl_lt(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 88 /* 1120 */:
            ShapeCoordinates::hexagon_bl_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 97 /* 1201 */:
            ShapeCoordinates::hexagon_tr_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 145 /* 2101 */:
            ShapeCoordinates::hexagon_lt_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 25 /* 0121 */:
            ShapeCoordinates::hexagon_lt_rb(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 70 /* 1012 */:
            ShapeCoordinates::hexagon_bl_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

        case 100 /* 1210 */:
            ShapeCoordinates::hexagon_bl_tr(cell, x0, x1, x2, x3, minV, maxV);
            break;

            /* 6-sided saddles */

        case 17 /* 0101 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::hexagon_lt_rb(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 68 /* 1010 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::hexagon_bl_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 153 /* 2121 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::hexagon_lt_rb(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 102 /* 1212 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::hexagon_bl_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

            /* 7-sided saddles */

        case 152 /* 2120 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 137 /* 2021 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_bl(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 98 /* 1202 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_tl(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 38 /* 0212 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 0 */
            if (center_avg == 2)
            {
                ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 18 /* 0102 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_tr(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 33 /* 0201 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_bl(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 72 /* 1020 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_tl(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 132 /* 2010 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            /* should never be center_avg === 2 */
            if (center_avg == 0)
            {
                ShapeCoordinates::triangle_br(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::heptagon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

            /* 8-sided saddles */

        case 136 /* 2020 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            if (center_avg == 0)
            {
                ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            else if (center_avg == 1)
            {
                ShapeCoordinates::octagon(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;

        case 34 /* 0202 */:
            center_avg = computeCenterAverage(x0, x1, x2, x3, minV, maxV);
            if (center_avg == 0)
            {
                ShapeCoordinates::tetragon_bl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_tr(cell, x0, x1, x2, x3, minV, maxV);
            }
            else if (center_avg == 1)
            {
                ShapeCoordinates::octagon(cell, x0, x1, x2, x3, minV, maxV);
            }
            else
            {
                ShapeCoordinates::tetragon_tl(cell, x0, x1, x2, x3, minV, maxV);
                ShapeCoordinates::tetragon_br(cell, x0, x1, x2, x3, minV, maxV);
            }
            break;
        }

        return cell;
    }

    void toGridData(const emscripten::val &jsData)
    {
        unsigned length = jsData["length"].as<unsigned>();
        std::cout << "length:" << length << "\n";

        std::vector<uint8_t> buffer(length);

        emscripten::val view(emscripten::typed_memory_view(
            buffer.size(),
            buffer.data()));

        view.call<void>("set", jsData);

        std::cout << (int)buffer[5] << "\n";
    }

    void test(const emscripten::val &jsData)
    {
        unsigned length = jsData["length"].as<unsigned>();
        std::cout << "length:" << length << "\n";

        std::vector<uint8_t> buffer(length);

        emscripten::val view(emscripten::typed_memory_view(
            buffer.size(),
            buffer.data()));

        view.call<void>("set", jsData);

        std::cout << (int)buffer[5] << "\n";
    }

    const emscripten::val isoBands(const emscripten::val &jsData,
                                   const int &width, const int &height,
                                   const emscripten::val &jsBands,
                                   const emscripten::val &jsGridBox,
                                   const emscripten::val &options)
    {
        bool flipY = false;
        bool flatten = false;
        int smoothIterations = 0;
        if (!options.isNull())
        {
            if (!options["flipY"].isNull())
            {
                flipY = options["flipY"].as<bool>();
            }
            if (!options["flatten"].isNull())
            {
                flatten = options["flatten"].as<bool>();
            }
            if (!options["smoothIterations"].isNull())
            {
                smoothIterations = options["smoothIterations"].as<int>();
            }
        }

        emscripten::val ret;

        std::vector<IsoBandResult> resultList;

        double west = jsGridBox[0].as<double>();
        double south = jsGridBox[1].as<double>();
        double east = jsGridBox[2].as<double>();
        double north = jsGridBox[3].as<double>();

        int bandNum = jsBands["length"].as<int>();

        std::vector<double> gridData;
        gridData.resize(height * width);

        emscripten::val view(emscripten::typed_memory_view(
            gridData.size(),
            gridData.data()));
        view.call<void>("set", jsData);

        std::vector<pair<double, double>> bandValues;

        std::vector<std::vector<std::vector<Cell>>> cellGridBreaks;

        for (int k = 0; k < bandNum; k++)
        {
            double minV = jsBands[k][0].as<double>();
            double maxV = jsBands[k][1].as<double>();
            bandValues.push_back({minV, maxV});

            std::vector<std::vector<Cell>>
                cellGrid;
            cellGrid.resize(width - 1);

            for (int i = 0; i < width - 1; ++i)
            {
                cellGrid[i].resize(height - 1);
                for (int j = 0; j < height - 1; ++j)
                {
                    cellGrid[i][j] = prepareCell(gridData, width, i, j, minV, maxV);
                }
            }

            cellGridBreaks.push_back(std::move(cellGrid));
        }

        for (int i = 0; i < cellGridBreaks.size(); i++)
        {
            auto [minV, maxV] = bandValues[i];
            IsoBandResult result;
            std::vector<IsoBandPolygon> polygons = traceBandPaths(gridData, width, height, cellGridBreaks[i], minV, maxV);

            std::sort(polygons.begin(), polygons.end(), [](const auto &a, const auto &b)
                      { return a.area > b.area; });

            std::vector<bool> allGroupBool;
            allGroupBool.resize(polygons.size());
            auto allGrouped = [&allGroupBool]()
            {
                for (size_t i = 0; i < allGroupBool.size(); i++)
                {
                    if (allGroupBool[i] == false)
                        return false;
                }
                return true;
            };

            while (!allGrouped())
            {
                for (size_t i = 0; i < polygons.size(); i++)
                {
                    if (!allGroupBool[i])
                    {
                        // create new group starting with the larger not already grouped ring
                        std::vector<IsoBandPolygon> group;
                        allGroupBool[i] = true;
                        group.push_back(std::move(polygons[i]));
                        // group all the rings contained by the outermost ring
                        for (size_t j = i + 1; j < polygons.size(); j++)
                        {
                            if (!allGroupBool[j])
                            {
                                if (isInside(polygons[j].coords, group[0]))
                                {
                                    allGroupBool[j] = true;
                                    group.push_back(std::move(polygons[j]));
                                }
                            }
                        }
                        // insert the new group
                        result.polygons.push_back(std::move(group));
                    }
                }
            }

            // std::cout<<"groupedLinearRingsCoords:"<<result.polygons.size()<<"\t"<<result.polygons[0].size()<<"\n";

            resultList.push_back(std::move(result));

            // std::cout << bandPolygons.size() << "\n";
            // ret.push_back(std::move(bandPolygons));
        }

        ret = emscripten::val::array();

        double originalWidth = east - west;
        double originalHeigth = north - south;
        double scaleX = originalWidth / (width - 1);
        double scaleY = originalHeigth / (height - 1);

        size_t polygonCount = 0;
        for (auto i = 0; i < resultList.size(); i++)
        {
            string minV = jsBands[i][0].call<emscripten::val>("toString").as<std::string>();
            string maxV = jsBands[i][1].call<emscripten::val>("toString").as<std::string>();

            auto properties = emscripten::val::object();
            properties.set("isoBand", minV + "-" + maxV);

            if (flatten)
            {
                const auto &bandPolygons = resultList[i].polygons;
                for (auto j = 0; j < bandPolygons.size(); j++)
                {
                    auto geometry = emscripten::val::object();
                    geometry.set("type", "Polygon");

                    auto ring = gen_polygon(bandPolygons[j], flipY,
                                          west, south, north, scaleX, scaleY, smoothIterations);
                    geometry.set("coordinates", ring);

                    auto feature = emscripten::val::object();
                    feature.set("type", "Feature");
                    feature.set("geometry", geometry);
                    feature.set("properties", properties);
                    ret.set(polygonCount++, feature);
                }
            }
            else
            {
                auto geometry = emscripten::val::object();
                geometry.set("type", "MultiPolygon");

                auto multiPolygon = emscripten::val::array();

                const auto &bandPolygons = resultList[i].polygons;
                for (auto j = 0; j < bandPolygons.size(); j++)
                {
                    auto ring = gen_polygon(bandPolygons[j], flipY,
                                          west, south, north, scaleX, scaleY, smoothIterations);
                    multiPolygon.set(j, ring);
                }

                geometry.set("coordinates", multiPolygon);

                auto feature = emscripten::val::object();
                feature.set("type", "Feature");
                feature.set("geometry", geometry);
                feature.set("properties", properties);

                ret.set(polygonCount++, feature);
            }
        }

        return ret;
    }
};
